﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.Statistics.RandomNumberGenerator
{
    /// <summary>
    /// Generates a random deviate from the beta distribution with parameters A and B. The density of the beta is 
    /// x^(a-1) * (1-x)^(b-1) / B(a,b).
    /// </summary>
    [Serializable]
    public class BetaRandomDeviateGenerator : AbstractRandomGenerator
    {
        /// <summary>
        /// The first parameter of the beta distribution.
        /// </summary>
        private readonly double a;

        /// <summary>
        /// The second parameter of the beta distribution.
        /// </summary>
        private readonly double b;

        /// <summary>
        /// This field holds the internal random number engine.
        /// </summary>
        private readonly Random engine;

        /// <summary>
        /// Initializes a new instance of the <see cref="BetaRandomDeviateGenerator"/> class.
        /// </summary>
        /// <param name="a">The first parameter of the beta distribution.</param>
        /// <param name="b">The second parameter of the beta distribution.</param>
        public BetaRandomDeviateGenerator(double a, double b)
            : base((uint) DateTime.Now.Millisecond)
        {
            if (a <= 0)
            {
                throw new ArithmeticException("a <= 0");
            }

            if (b <= 0)
            {
                throw new ArithmeticException("b <= 0");
            }

            this.a = a;
            this.b = b;
            this.engine = new Random((int) this.Seed);
        }

        /// <summary>
        /// Initializes a new instance of the <see cref="BetaRandomDeviateGenerator"/> class.
        /// </summary>
        /// <param name="a">The first parameter of the beta distribution.</param>
        /// <param name="b">The second parameter of the beta distribution.</param>
        /// <param name="seed">A number used to calculate a starting value for the pseudo-random number sequence.
        /// If a negative number is specified, the absolute value of the number is used.</param>
        public BetaRandomDeviateGenerator(double a, double b, uint seed)
            : base(seed)
        {
            if (a <= 0)
            {
                throw new ArithmeticException("a <= 0");
            }

            if (b <= 0)
            {
                throw new ArithmeticException("b <= 0");
            }

            this.a = a;
            this.b = b;
            this.engine = new Random((int) seed);
        }

        /// <summary>
        /// Returns a nonnegative random number.
        /// </summary>
        /// <returns>
        /// A 32-bit signed integer greater than or equal to zero and less than maximum integer value.
        /// </returns>
        public override int Next()
        {
            return this.Next(0, int.MaxValue);
        }

        /// <summary>
        /// Nexts the specified max value.
        /// </summary>
        /// <param name="maxValue">Returns a nonnegative random number less than the specified maximum.</param>
        /// <returns>
        /// A 32-bit signed integer greater than or equal to zero and less than maximum integer value.
        /// </returns>
        public override int Next(int maxValue)
        {
            return this.Next(0, maxValue);
        }

        /// <summary>
        /// Returns a random number within a specified range.
        /// </summary>
        /// <param name="minValue">The inclusive lower bound of the random number returned.</param>
        /// <param name="maxValue">The exclusive upper bound of the random number returned. maxValue must be
        /// greater than or equal to minValue.</param>
        /// <returns>
        /// A 32-bit signed integer greater than or equal to zero and less than maximum integer value.
        /// </returns>
        public override int Next(int minValue, int maxValue)
        {
            if (minValue < 0)
            {
                minValue = 0;
            }

            if (maxValue == int.MaxValue)
            {
                maxValue -= 1;
            }

            if (minValue > maxValue)
            {
                int tempuri = maxValue;

                maxValue = minValue;
                minValue = tempuri;
            }

            double number = this.GenerateRandomNumber(this.a, this.b);

            if (number < 0)
            {
                number = Math.Abs(number);
            }

            if (number > 1)
            {
                number = number - Math.Truncate(number);
            }

            return (int) (Math.Floor((maxValue - minValue + 1)*number + minValue));
        }

        /// <summary>
        /// Returns a random number between 0.0 and 1.0.
        /// </summary>
        /// <returns>
        /// A double-precision floating point number greater than or equal to 0.0, and less
        /// than 1.0.
        /// </returns>
        public override double NextDouble()
        {
            return this.GenerateRandomNumber(this.a, this.b);
        }

        private double GenerateRandomNumber(double aa, double bb)
        {
            const double expmax = 89.0;
            const double infnty = 1.0e38;
            double paramA = 0;
            double paramB = 0;
            double olda = -1.0;
            double oldb = -1.0;
            double alpha = 0;
            double beta = 0;
            double gamma = 0;
            double k1 = 0;
            double k2 = 0;
            double r, s, t, u1, u2, v, w, y, z;
            double genbet;

            bool qsame = ((olda == aa) && (oldb == bb));

            if (qsame)
            {
                goto S20;
            }

            if (!(aa <= 0.0 || bb <= 0.0))
            {
                goto S10;
            }

            S10:
            olda = aa;
            oldb = bb;
            S20:
            if (!(Math.Min(aa, bb) > 1.0))
            {
                goto S100;
            }

            if (qsame)
            {
                goto S30;
            }

            paramA = Math.Min(aa, bb);
            paramB = Math.Max(aa, bb);
            alpha = paramA + paramB;
            beta = Math.Sqrt((alpha - 2.0)/(2.0*paramA*paramB - alpha));
            gamma = paramA + 1.0/beta;
            S30:
            S40:
            u1 = this.engine.NextDouble();
            u2 = this.engine.NextDouble();
            v = beta*Math.Log(u1/(1.0 - u1));

            if (!(v > expmax))
            {
                goto S50;
            }

            w = infnty;
            goto S60;
            S50:
            w = paramA*Math.Exp(v);
            S60:
            z = Math.Pow(u1, 2.0)*u2;
            r = gamma*v - 1.3862944;
            s = paramA + r - w;

            if (s + 2.609438 >= 5.0*z)
            {
                goto S70;
            }

            t = Math.Log(z);

            if (s > t)
            {
                goto S70;
            }

            if (r + alpha*Math.Log(alpha/(paramB + w)) < t)
            {
                goto S40;
            }

            S70:

            if (!(aa == paramA))
            {
                goto S80;
            }

            genbet = w/(paramB + w);
            goto S90;
            S80:
            genbet = paramB/(paramB + w);
            S90:
            goto S230;
            S100:

            if (qsame)
            {
                goto S110;
            }

            paramA = Math.Max(aa, bb);
            paramB = Math.Min(aa, bb);
            alpha = paramA + paramB;
            beta = 1.0/paramB;
            double delta = 1.0 + paramA - paramB;
            k1 = delta*(1.38889E-2 + 4.16667E-2*paramB)/(paramA*beta - 0.777778);
            k2 = 0.25 + (0.5 + 0.25/delta)*paramB;
            S110:
            S120:
            u1 = this.engine.NextDouble();
            u2 = this.engine.NextDouble();

            if (u1 >= 0.5)
            {
                goto S130;
            }

            y = u1*u2;
            z = u1*y;

            if (0.25*u2 + z - y >= k1)
            {
                goto S120;
            }

            goto S170;
            S130:
            z = Math.Pow(u1, 2.0)*u2;

            if (!(z <= 0.25))
            {
                goto S160;
            }

            v = beta*Math.Log(u1/(1.0 - u1));

            if (!(v > expmax))
            {
                goto S140;
            }

            w = infnty;
            goto S150;
            S140:
            w = paramA*Math.Exp(v);
            S150:
            goto S200;
            S160:

            if (z >= k2)
            {
                goto S120;
            }
            S170:
            v = beta*Math.Log(u1/(1.0 - u1));

            if (!(v > expmax))
            {
                goto S180;
            }

            w = infnty;
            goto S190;
            S180:
            w = paramA*Math.Exp(v);
            S190:
            if (alpha*(Math.Log(alpha/(paramB + w)) + v) - 1.3862944 < Math.Log(z))
            {
                goto S120;
            }

            S200:
            if (!(paramA == aa))
            {
                goto S210;
            }

            genbet = w/(paramB + w);
            goto S220;
            S210:
            genbet = paramB/(paramB + w);
            S230:
            S220:

            return genbet;
        }
    }
}